Long-time behavior of the momentum distribution during the sudden expansion of 
a spin-imbalanced Fermi gas in one dimension 
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We study the sudden expansion of spin-imbalanced ultracold lattice fermions with attractive interactions in 
one dimension after turning off the longitudinal confining potential. We show that the momentum distribution 
functions of majority and minority fermions quickly approach stationary values due to a quantum distillation 
mechanism that results in a spatial separation of pairs and majority fermions. As a consequence, Fulde-Ferrell- 
Larkin-Ovchinnikov (FFLO) correlations are lost during the expansion. Furthermore, we argue that the shape of 
the stationary momentum distribution functions can be understood by relating them to the integrals of motion in 
this integrable quantum system. We discuss our results in the context of proposals to observe FFLO correlations, 
related to recent experiments by Liao et ai, Nature 467, 567 (2010). 
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The combination of strong correlations and quantum fluctu- 
ations makes one-dimensional (ID) systems the host of exotic 
phases and physical phenomena llj,|2|]. Those phases and phe- 
nomena, in many occasions first predicted theoretically, have 
been observed in condensed matter experiments and have be- 
gun to be studied with ultracold atomic gases J^. A sys- 
tem of particular interest in recent years has been the spin 
imbalanced ID Fermi gas. Following theoretical predictions 
(JO], its grand canonical phase diagram has recently been 
investigated experimentally ifioll . The major interest in this 
model comes from the fact that its entire partially polarized 
phase has been theoretically shown Iblld ll ll4l5ll (for a review, 



is to let the gas expand in the ID lattice after turning off the 
longitudinal confining potential, and then measure the den- 
sity profiles or the MDFs of the independent species and/or 
pairs after some expansion time. Some aspects of such an 
expansion experiment have already been successfully carried 



out in ID tubes 0(1 13111 as well as in 2D and 3D optical lat- 



see JljJ) t° be the ID-analogue of the Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) state 1171 1 1811 . The FFLO phase was 
introduced to describe a possible equilibrium state in which 
magnetism and superconductivity coexist thanks to the forma- 
tion of pairs with finite center-of-mass momentum leading to a 
spatially oscillating order parameter. The existence of such a 
phase has remained controversial in dimensions higher than 
one in theoretical studies 1 1^-21], while experiments have 
found no evidence of the FFLO phase in three-dimensional 
systems l22l 12311 . 

An important challenge in ultracold fermion experiments, 
which may have already realized the FFLO state lHoll . is to 
confirm the existence of FFLO correlations (for recent pro- 
posals see, e.g., J24-27]). A direct measurement of the pair 
momentum distribution function (MDF) in the partially po- 
larized state [3 has been suggested to provide such 
evidence J28ll . However, this remains very difficult because 
after turning off all confining potentials, the transverse expan- 
sion (in the directions of very tight confinement) dominates 
over the longitudinal one j29l]. Another interesting possibility 



tices 0211 . namely the independent control over lattice and the 
trapping potential and the measurement of the density pro- 
files after the expansion. For ID gases, interaction effects 
during the expansion cannot in general be neglected, leading 
to fundamentally different behavior of observables before and 
after the gas has expanded. For example, the expansion of 
the Tonks-Girardeau gas in ID results in a bosonic gas with 



a fermionic MDF 133 



g) states of bosons 1 3a 



35 1, and initially incoherent (insulat 



3711 and fermions 113 811 can develop 



quasi-long range correlations during the expansion. 

The question we are set to address is the fate of the MDFs of 
fermions and pairs during an expansion in one dimension, as 
described by the attractive Hubbard model. We use a combi- 
nation of numerical simulations, based on the time-dependent 
density matrix renormalization group approach (t-DMRG) 
[39l|40|], and analytical (Bethe-Ansatz) results. We first show 
that the MDFs of majority and minority fermions become sta- 
tionary after a relatively short expansion time, t ~ La/ J, 
where Lq is the initial size of the cloud and J is the hopping 
amplitude. For strong interactions, we explain this behavior in 
terms of a quantum distillation process 114 IB . as a consequence 
of which FFLO correlations are destroyed during the expan- 
sion. Finally, we discuss how these stationary MDFs can be 
theoretically understood within the framework of the Bethe- 
Ansatz. Our results suggest that the final form of the MDFs of 
minority and majority fermions are related to the distributions 



2 



of Bethe-Ansatz rapidities (a full set of conserved quantities) 
of this integrable lattice system. 



0.04 r 



U=-10J, N=8, p=0.5 W 



The Hubbard model (in standard notation 114210 reads: 



L-l 



H = -jJ2 (4+i,^,c + H.c.) + uj2 



(1) 



As the initial state, we always take the ground state of a 
trapped system. In the main text, we focus on a box trap, 
i.e., particles confined into a region of length Lq while we 
present results for the expansion from a harmonic trap in the 
supplementary material 14311 . We study lattices with L sites, 
N particles, and a global polarization of p = (N^ — N±)/N, 
where N a = J2e( n (°-)- All positions are given in units of 
the lattice spacing and momenta in inverse units of the lattice 
spacing (h = 1). 

The expansion is triggered by suddenly turning off the con- 
fining potential, thus allowing particles to expand in the lat- 
tice. We then follow the time-evolution using the numerically 
exact f-DMRG algorithm JHH- 

We use a Krylov-space 
based time-evolution method and enforce discarded weights 
of 10~ 4 or smaller with a time-step of St = 0.25/ J. Our 
main focus is on the time-evolution of the three MDFs: the 
ones for majority (er = f) and minority fermions (er = J,), de- 
noted by rik,cr and the pair MDF, Uk,p- These functions are 
computed from the corresponding one-particle (A =t, 4) or 
one-pair (A = p) density matrices via a Fourier transform 



(2) 



— w t ^ e j_ and A stands for j, 4-, P- We 



where V] )(T = c] j(7 ,^| 

normalize the MDFs so that J2 k Uk,\ = N\ (note that N p — 
^ £ (n^n^), i.e., it is equal to the total double occupancy in 
the system). 

For the expansion from a box, we focus on an initial density 
fixed ton = N j Lq = 0.8. In our f-DMRG simulations, which 
were carried out for N — 8 and N = 16 (Lq = 10 and 20, 
respectively) and various values of U, we were able to reach 
times of order t max ~ 80/ J for large U and t max ~ 40/ J for 
intermediate values of U ~ 4 J. t max also depends on p, with 
small values of p being more demanding. 

Typical results for the three MDFs of interest are presented 
in Fig. Q]for U = —10 J and p = 0.5 (corresponding to 
JVf = 6 and N± = 2; see the supplementary material for more 
results 14311 ). During the time evolution, they are all seen to 
quickly approach time-independent forms. In Fig. |TJa), it is 
apparent that the MDF of the majority fermions becomes nar- 
rower and develops small oscillations in the vicinity of k = 
as time passes. We find that those oscillations become smaller 
in amplitude and get restricted to smaller values of k after long 
expansion times, i.e., they seem to be a transient feature not 
present in the asymptotic distributions. The momentum distri- 
bution of the minority fermions [Fig.QJb)], on the other hand, 
becomes broader during the time evolution. 

The time evolution of the MDF of the pairs, depicted in 
Fig. Qlc), yields information on the fate of FFLO correlations 
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FIG. 1: (Color online) MDF for the expansion from a box trap (U = 
— 10 J, N — 8, p — 0.5, Lo — 10): (a) Uk,^, (b) nk,i, and (c) nk, p . 
The insets show the difference A\ (A =f , 4-, p, see text) between the 
MDF at a time t compared to the one at the largest time reached in 
the simulation. The vertical lines in the main panel in (c) mark the 
position of the FFLO wave-vector Q = ±nnp. 



in the expanding cloud. In the FFLO state, nu, P has maxima 
at Q = ±(kF j t — kFi) H]. These are visible in the t = curve 
(dashed line), where ±Q are marked by vertical lines. As the 
comparison of nk, p (t > 0) with the initial nk :P {t = 0) shows, 
the peaks at ±Q rapidly disappear, and nk. P (t) becomes nar- 
rower. In addition, new and shallower peaks form at fc < Q. 
Since we do not find those peaks at the same values of fc for 
other values of when N /Lo and p are the same, and we 
do not find them for all values of U, N/ Lq, and p studied, 
they appear to be related to finite-size effects. Hence, the dou- 
ble peak structure in nk, P (t — 0), which makes evident the 
presence of FFLO correlations in the initial state, is found to 
disappear during the expansion. Even though the FFLO corre- 
lations are lost during the expansion, the integral over the pair 
MDF, which equals the total double occupancy, does not van- 
ish. This implies that not all interaction energy is converted 
into kinetic energy and that some fraction of the original pairs 
remains by the time the MDFs have become stationary, which 
in experiments could be probed by measuring the double oc- 
cupancy. 

In order to quantify how the three MDFs above approach 
stationary forms, in the insets in Fig. [T] we plot A.\(t) = 
J2k \ n k,\(t) - nk,x(tm a x)\/Y,k n k,x(tma,x) vs t. These re- 
sults make apparent that the approach is close to exponen- 
tial for rikrf and rik,\. [insets in Fig. [Ha) and [Tib)], while it 
is power law for rik, P [inset in Fig.[TJc)] HH]. Remarkably, 
for the parameters of Fig. [TJ already at t J ~ 10, all Aa are 
< 10%. This means that the stationary MDFs obtained in 
this work should be achievable in current optical lattice setups 
1 3211 . A comparison between expansions from different box 
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as opposed to exponential, relaxation observed for nk, p (t) in 
Fig.Ec). 
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FIG. 2: (Color online) Natural orbital |3>o| corresponding to the 
largest eigenvalue of the pair-pair correlator P(£,j) (dashed lines) 
and spin density (SI ) (solid lines), (a) t — 0, (b) tj = 10. These 
results are for U = -10J, L = 20, N = 16 and p = 0.75, 
corresponding to Nf — 14 and N± = 2. 



sizes suggests that the emerging time scale in the observables 
with exponential relaxation is proportional to Lq. The origin 
of that time scale will be discussed below. 

While we are focusing the discussion on the case of the 
expansion from a box trap, we stress that the results for the 
MDFs in an expansion from a harmonic trap are quite sim- 
ilar (for an example, see the supplementary material 043J). 
Namely, we observe a comparably fast convergence of the 
MDFs to a stationary form and the disappearance of the peaks 
at ±Q in n^ p . The latter indicates the disappearance of FFLO 
correlations. 

To understand how the FFLO state breaks down as the gas 
expands, we calculate the eigenvector <&o of the pair-pair cor- 
relator P(£, m) = (ipg p ip m ,p) that corresponds to the largest 
eigenvalue. |$ |, shown in Fig. Eta), unveils the spatial struc- 
ture of the quasi-condensate in the initial state: it has an oscil- 
latory behavior with nodes (see also Ref. J5|). In these nodes, 
the spin density has its maxima to accommodate the major- 
ity fermions (Fig. |2t a ); see also 14311 ). indicative of the spin- 
density wave character with a modulation of (2Q)" 1 in the 
FFLO state. During the expansion, the nodes in |<5rj| disappear 
while |$q| develops a maximum at L/2, exceeding its initial 
value [see Fig. Efb)] . The latter is a consequence of a quan- 
tum distillation mechanism, described in Ref. 14 ill for U > 0, 
which allows the unpaired fermions to move away from the 
center of the system (i.e., they escape from the nodes of 
\3>o(t = 0)|). Loosely speaking, during first-order processes 
unpaired fermions exchange their positions with the pairs (a 
minority fermion hops towards the center of the trap), allow- 
ing the former to expand while the pairs move towards the 
center of the trap. This occurs over a time scale proportional 
to Lq and inversely proportional to J, which explains the time 
scale observed in the exponential approach of the majority 
and minority fermions to their stationary values. Once the 
unpaired fermions have spatially separated themselves from 
the pairs, they form a non-interacting gas whose MDF is sta- 
tionary. On much longer time scales (assuming U > 4 J), 
we expect the pairs to slowly expand as well. This transient 
dynamics of the pairs may be the reason for the power-law, 



In a recent work |45|], extrema in the spin-density of the 
expanding gas were observed in numerical calculations using 
various approaches. By comparing with the time-dependence 
of the order parameter within a time-dependent Bogoliubov- 
deGennes approach, it was argued that they are related to 
FFLO correlations. Our results show that, in a lattice sys- 
tem, the nodal structure of the FFLO state is ultimately lost 
as the system expands. Note, however, that in Ref. 14511 the 
main focus was on rather small polarizations p 13 ^ lead- 
ing to a wide partially polarized core before the expansion. 
We therefore expect the quantum distillation mechanism to 
take much longer to depolarize the core than what has so far 
been reached in numerical simulations |45|], leaving this case 
as an open question. 

We are now in a position to explain the antic orrelated be- 
havior of rik,f and mentioned in the discussion of Fig.Q] 
For large values of U, N p is essentially equal to iVj. and is ap- 
proximately unchanged during the expansion, rendering the 
interaction energy almost time independent. This implies that 
also the kinetic energy E^ m = — 2 JJ2k cos ^( n fc,t + n k.f) 
is approximately conserved, which is only possible if the two 
MDFs behave in the opposite way during the expansion. The 
broadening of the minority MDF rik.i with respect to the ini- 
tial state is a direct consequence of the spatial separation of 
excess fermions from the pairs, leaving the latter confined in 
the center of the cloud. Since in the center the local polar- 
ization decreases, the stationary form of nk,\. is well approxi- 
mated by the equilibrium one for equal populations Nf = N± 
instead of Nf > N l H. 

The fact that the MDFs become stationary after the expan- 
sion from a box or a harmonic trap is in itself not surprising, as 
in the limit of long expansion times, the cloud becomes very 
dilute with, for the attractive case, the typical inter-particle 
distance being much larger than the bound-state size. Hence, 
one may assume that pairs and unpaired particles are essen- 
tially noninteracting. The MDF in such an asymptotic limit 
should be determined by the initial conditions right after the 
quench. For instance, for generic models, the total energy 
(which is conserved during the expansion) plays a fundamen- 
tal role in determining the expansion dynamics (see Ref. ll46tl 
for a related work for U > 0). For an integrable model, such 
as the (attractive) Hubbard model of Eq. Q}, all integrals of 
motion are in principle known from the Bethe Ansatz and are 
conserved during the expansion |42J]. We argue below how 
to interpret the shape of certain stationary MDFs in terms of 
such integrals of motion. This is closely related to the previ- 
ously studied fermionization of the MDF of an expanding gas 
of hard-core bosons 03K35 1 . 

For the model studied here, we first note that the for- 
mation of a distinct minimum in the difference distribution 
Silk = nk,f — nk.i [see Figs. [3ja) and (b)] is reminiscent of 
the corresponding distribution of real-valued charge rapidities 
(for intermediate U) in the ground state in a box. From the 
point of view of the rapidity distributions, they need to be de- 



4 



termined right after turning off the trap and the subsequent 
expansion does not play any role; it is the MDFs which will 
evolve and asymptotically approach the former as the expan- 
sion proceeds 14711 . We can calculate the pre-quench values of 
the rapidities by numerically solving the Bethe-Ansatz equa- 
tions for a system of size Lq and open boundary conditions 



114 8145011 . For the ground state of the attractive Hubbard model, 
there are two types of rapidities present: real- and complex- 
valued charge rapidities («„ and K a ) which correspond to un- 
paired fermions and pairs, respectively (y = 1, . . . , — iVj,, 



1. 



, 2AT|_, with n a and n* a appearing pairwise). 



To calculate the effect of the quench of the trapping poten- 
tial exactly is in principle possible but complicated in prac- 
tice 115 ill , so we will resort to some simplifications. To start, 
we assume that the number of pairs is conserved during the 
quench, and thus no pure-spin excitations are produced. Fur- 
ther, we use the observation that the overlap between the pre- 
quench eigenstate and the post-quench state has a maximum 
amplitude for components of the latter with the same set of 
rapidities pill . We then identify, asymptotically, the distri- 
bution of real-valued charge rapidities with that of unpaired 
fermions (Srik), and of the real part of complex-valued (string) 
charge rapidities with that of minority fermions (n^.i) since 
they remain paired. Finally, we model the quench by convolv- 
ing the pre-quench distributions p\ = (1/2) 6(k ± K v ) 
and p2 = (1/2) ^ R.ck ct ) with the (periodized) ker- 

nels: (i) Lq sinc 2 (fcio/2) for the former and (ii) a simple 
Lorentzian for the latter. The first choice is inspired by the 
exact result for the release of a single particle from a box, 
while the second choice is done for simplicity given that the 
results are relatively featureless in comparison. Illustrative re- 
sults are shown in Fig. [3] and the agreement is very good, spe- 
cially away from the Brillouin-zone center. Note that there 
are no fitting parameters in the case of 5rik and a single fitting 
parameter, the width of the Lorentzian, in the case of nk,i- 

In conclusion, we demonstrated that the initial FFLO state 
is destroyed during the expansion of an attractively interacting 
partially polarized ID Fermi gas, and that direct signatures 
of the FFLO phase in the initial pair MDF are washed out 
as a consequence of interactions. Nevertheless, the sudden 
expansion is an interesting non-equilibrium experiment that 
through the asymptotic form of the MDFs yields information 
on the initial state. Our analysis suggests that the shape of the 
MDFs can be related to the distribution of rapidities, which 
constitute a full set of integrals of motion for this integrable 
quantum model and fully determine the initial state. Since we 
showed that the MDFs of majority and minority fermions as 
well as the one of pairs rapidly take a stationary form, this 
should be accessible on typical experimental time-scales. 
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FIG. 3: (Color online) Comparison of the stationary MDFs Snk = 
Tlk,t — n k,i [(a),(b)] and nk,i [(c), (d)] for the expansion from a box 
with N = 8, p = 0.5 (corresponding to N t = 6, N± = 2) [(a),(c): 
U = — 4J, (b),(d): U = — 10 J] to the form expected from the ra- 
pidities known from the Bethe-Ansatz: i-DMRG (solid lines), mod- 
els discussed in the text (dashed lines). The vertical lines mark the 
positions of the rapidities. 
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Additional results for the MDFs of a spin-imbalanced Fermi gas 
with attractive interactions 



We here provide additional t-DMRG results for the MDFs 
nk,x °f a spin-imbalance Fermi gas with attractive interac- 
tions, expanding from a box trap. Figure [4] contains data for 
N = 16 (with U = -10J, p = 0.75, L = 20). In this 
case we were able to reach maximum times of t max ~ 30/ J. 
Nevertheless, a fast approach towards a stationary form is ob- 
vious from this figure, corroborating the conclusions of the 
main text (see the discussion of Fig. 1 of the main text). The 
same applies to the qualitative trends: rik,f shrinks while rik,\. 
broadens. 

In Fig. [5] we display results for U = — 4 J and N = 8 with 
a polarization of p = 0.5. In this case, the convergence to 
a stationary form is evident in all three MDFs. Note that, in 
contrast to the case of U = —10 J discussed in the main text, 
nfe.t does not exhibit any transient fluctuations at small fc. 




FIG. 4: MDFs for (a) spin up, (b) spin down, and (c) pairs for the 
expansion from a box trap with U = —10 J, N = 16, p — 0.75, 
L = 20, plotted at times tj = 0, 10, 20, 30. The circles in the 
central panel represent the MDF of an unpolarized gas with Nf = 

Ni = 2. 



0.04 



0.02 - 



i 0.01 4 



0.01- 




FIG. 5: MDFs for (a) spin up, (b) spin down, and (c) pairs for 
the expansion from a box trap with U = — 4J, N = 8, p = 0.5, 
I/O = 10, plotted at times t J = 0, 20, 40. The circles in the central 
panel represent the MDF of an unpolarized gas with Nf = N± = 2. 



separate from unbound fermions, such that the cloud develops 
a two-shell structure with a fully paired core and fully polar- 
ized wings containing the excess Nf — N± fermions. This sug- 
gests that the asymptotic momentum distribution of the mi- 
nority component should be approximated by its ground-state 
value before the expansion calculated in the absence of excess 
fermions, that is for Nf — N±. The results for the MDF rik,i 
that we obtain using this assumption are plotted in Figs.|4]and 
[5] with circles. Indeed, we see a rather good agreement with 
the stationary form of the MDF n^, where the latter was 
calculated with i-DMRG. In particular, in the limit of large 
initial polarization p — > 1, the number of pairs is very small. 
In this low density (or equivalently, strong-coupling) regime 
the ground state momentum distribution becomes equal to 
the square of the Fourier transform of the molecular wave- 
function for the relative motion: 



Discussion of the qualitative behavior of the MDFs of the 
spin-imbalanced Fermi gas with attractive interactions 



\u\ 3 L 

VU 2 + 16J 2 (-4 J cos k + yJU 2 + 16 J 2 ) 2 



(3) 



Comparing Figs. l(a)-(b) of the main text, as well as 
Fig. |3 a) and (b) shown here in the supplementary material, 
we see that the momentum distribution n^.f of the majority 
component shrinks during the expansion, whereas the distri- 
bution nk,i of the minority component broadens significantly. 
This is the result of collisions between up and down fermions 
which take place in the inner part of the system, where both 
spin components are present, and transfer momenta between 
them. In the long time limit and for \ U\ > J, the pairs phase 



The corresponding shrinking of the majority momentum dis- 
tribution during the expansion can then be understood from 
conservation of total energy. Indeed, for | U\ ^> J the num- 
ber of double occupancies remains close to N p = N± dur- 
ing the expansion, implying that the interaction energy in the 
Hubbard model is essentially time independent. As a conse- 
quence, the kinetic energy E^ n = — 2 J^2k cos ^( ri fe,t+ ri fc4) 
is also conserved, implying that the distribution rik.f must 
shrink to compensate the broadening of rik.i- 
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FIG. 6: Natural orbital |<E>o| corresponding to the largest eigenvalue 
of the pair-pair correlator P(£,j) (dashed lines) and spin density 
(Sf) (solid lines). (a),(c) t = 0, (b),(d) tj = 10. These results 
are for U = -10 J, L = 20, N = 16 and p = 0.25 [panels (a) and 
(b)], and p=0.5 [panels (c) and (d)]. 



Additional results for the time-evolution of the quasicondensate 

In Fig. 2 of the main text, we show the spin-density (Sf) 
and the eigenvector |<f>o| of the pair-pair correlator in the 
initial state and after an expansion time of t J = 10 for 
U = —10 J and p = 0.75 (corresponding to Nf = 14 and 
JVj. = 2 with Lq = 20). Here we present additional results for 
p = 0.25 (JV t = 10 and N± = 6) and (JV t = 12 and N t = 4) 
in Fig. [6] The nodal structure of quasicondensate is best seen 
in the case of p = 0.25 [Fig.|6|a)] with the spin-density taking 
maxima in the nodes of |$n|. During the expansion, we ob- 
serve the same spatial demixing of excess fermions from the 
pairs that is discussed in the main text for the case ofp = 0.75. 
Note that the demixing occurs over short expansion times dur- 
ing which the cloud has expanded by about a factor of two 
only. 

In Fig. [7] we demonstrate that the demixing is not limited 
to the strongly interacting regime. We present results for (Sf ) 
and |$ | for U = -4 J at p = 0.75 with N = 8. Obvi- 




FIG. 8: (Color online) MDFs for the expansion from a harmonic 
trap: (a) rik,t, (b) nk,i and (c) n^p.These results were obtained for 
N = 8, U = -10J, p = 0.75, V = 0.016J, and at times tj = 
0,90, 100. Inset in (a): initial density (m) (solid line) and spin- 
density profile 2{Sf) (dashed line). 



ously, at time t J = 10, we observe that the spatial structure 
of the FFLO-quasicondensate is lost and that excess fermions 
and pairs are separated from each other. The main reason for 
the separation of excess fermions from pairs is, as discussed 
in the main text, the difference in the bare velocities of pairs 
and excess fermions. In the lattice this difference is large at 
U ^> J yet the quantum distillation mechanism still works 
even at intermediate values of \U\ ~ 4 J as shown here. The 
only noticeable difference between large U and intermediate 
U is that the density of the pairs in the core of the cloud does 
not increase as the system expands at U = — 4J. Based on 
our results for U = —4 J, we, therefore, speculate that in the 
contiuuum, where the bare velocities of excess fermions ver- 
sus pairs differ by a factor of two only, the demixing of excess 
fermions and pairs should also occur during the expansion. 
However, it is numerically very demanding to simulate this 
limit using time-dependent DMRG since very low densities 
and therefore, long expansion times would be necessary. 
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FIG. 7: Natural orbital | <E>o | corresponding to the largest eigenvalue 
of the pair-pair correlator P(l,j) (dashed lines) and spin density 
(Sf ) (solid lines), (a) t = 0, (b) t J = 10. These results are for 
U = -4J, L = 10, N = 8 and p = 0.75. 



MDFs for the sudden expansion from a harmonic trap 

In relation with experiments, it is also important to incorpo- 
rate the harmonic confinement, i.e., H trav = Vq X)/=i n e 
L/2) 2 . To that end, we have prepared a spin-imbalanced sys- 
tem with U = — 10 J in a harmonic trap with Vq > Ofori < 0, 
and then quenched the trapping potential to Vq = at t = 0. 
For the parameters of Fig. [8] the partially polarized phase that 
sits in the core is surrounded by fully polarized wings (see the 
inset in Fig. [8). During the expansion, one can see that the 
behavior of the MDFs is very similar to the one starting from 
a box in Fig. 1 of the main text. All MDFs become stationary 
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FIG. 9: MDF nu, a = n,k,t = nu.i for the expansion from a box 
trap with U = 8 J, N = 8, p = 0, plotted at times tJ = 0, 5, 10, 20. 
We observe that in the case of repulsive interactions, much shorter 
times can be reached than for U < 0. On the accessible time scales, 
the MDF still changes, yet at its edge, the curves for t J = 10 and 
t J = 20 lie on top of each other. 



shortly after the release from the trap. The stationary nj, -j- is 
narrower while rik.i is broader than their corresponding initial 
distributions, and the double peak structure in Uk, v , which is 
due to the FFLO correlations in the initial state, disappears. 



symmetries as well, restricting the analysis to p = 0, allowing 
us to reach t ~ 25/ J for U = 8 J (see Fig. [9). In the case of 
U = 8 J, there are no pairs, and hence over the full extent of 
the expanding cloud, majority and minority fermions can still 
interact, whereas in the case of U < and p > 0, fast ma- 
jority fermions escape [3] and the pairs and majority fermions 
mostly decouple due to the quantum distillation mechanism 
that is described in the main text. This is likely the reason 
why in the repulsive gas with U > and p = 0, there is a 
stronger entanglement during the expansion. On the accessi- 
ble time scales, the MDF nfc 1<7 of the repulsive gas still un- 
dergoes changes, yet in the edge of the MDF, the curves at 
the longest times coincide (see also (H). The case of U > 
thus sets an example where the quantum simulation with ultra- 
cold atomic gases could help us to go to longer times than 
what is currently possible with numerical methods to clarify 
the asymptotic behavior of the MDF (compare the relaxation 
dynamics problem studied in Ref. J3l). Note that for the ex- 
pansion of a repulsive gas with initial densities (n*) < 1, the 
double-occupancy decreases in contrast to the attractive 
case, discussed in the main text. In the attractive case, the sur- 
vival of a certain fraction of the initial double occupancy is 
expected due to the presence of pairs. 



Time-evolution of the MDFs of a two-component Fermi gas with 
repulsive interactions 

We have also studied the time-evolution of other ID mod- 
els during the expansion, including most notably the repulsive 
Hubbard model with p = (compare HI). 

The U > case turns out to be a numerically much harder 
problem for t-DMRG, as entanglement grows much faster 
(see the review J2l for how entanglement growth limits t- 
DMRG). Therefore, we resorted to exploiting non-Abelian 



[1] F. Heidrich-Meisner, M. Rigol, A. Muramatsu, A. E. Feiguin, 

and E. Dagotto, Phys. Rev. A 78, 013620 (2008). 
[2] U. Schollwock, Ann. Phys. (NY) 326, 96 (2011). 
[3] J. Kajala, F. Massel, and P. Torma, Phys. Rev. A 84, 041601(R) 

(2011) . 

[4] S. Trotzky, Y.-A. Chen, A. Flesch, I. P. McCulloch, 
U. Schollwock, J. Eisert, and I. Bloch, Nature Phys. 8, 325 

(2012) . 



